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On the probability that all eigenvalues of 
Gaussian and Wishart random matrices lie 

within an interval 

Marco Chiani, Fellow, IEEE 

Abstract 

We derive the probability 'tp{a,b) = Pr(a < Amin(M), Aniax(M) < b) that all eigenvalues of a 
random matrix M lie within an arbitrary interval [a, b], when M is a real or complex finite dimensional 
Wishart, double Wishart, or Gaussian symmetric/hermitian matrix. We give efficient recursive formulas 
allowing, for instance, the exact evaluation of ^/;(a, b) for Wishart matrices with number of variates 
500 and degrees of freedom 1000. We also prove that the probability that all eigenvalues are within 
the limiting spectral support (given by the Marchenko-Pastur or the semicircle laws) tends for large 
dimensions to the universal values 0.6921 and 0.9397 for the real and complex cases, respectively. 
Applications include improved bounds for the probability that a Gaussian measurement matrix has a 
given restricted isometry constant in compressed sensing. 

Index Terms 

Random Matrix Theory, Principal Component Analysis, Compressed Sensing, eigenvalues distribu¬ 
tion, Tracy-Widom distribution, Wishart matrices, Gaussian Orthogonal Ensemble, MANOVA. 

I. Introduction 

Many mathematical models in information theory and physics are formulated by using matrices 
with random elements. In particular, the distribution of the eigenvalues of Wishart and Gaussian 
random matrices plays a key role in multivariate analysis, including principal component analysis, 

Marco Chiani is with DEI, University of Bologna, V.le Risorgimento 2, 40136 Bologna, ITALY (e-mail: 
marco.chiani@uniho.it). This work was partially supported hy the Italian Ministry of Education, Universities and Research 
(MIUR) under the PRIN Research Project “GRETA”. 


Submitted 2015 


DRAFT 



VOL. X, NO. Y, MONTH 2015 


2 


analysis of large data sets, information theory, signal processing and mathematical physics 
[1]-[11]. Most of the problems concern the distribution of the smallest and/or of the largest 
eigenvalue, or of a randomly picked (unordered) eigenvalue. 

For example, in compressed sensing the probability that a randomly generated measurement ma¬ 
trix has a given restricted isometry constant is related to the probability Pr(a < Amin(M), Amax(M) < 
b), where Amin(M), Amax(M) denote the minimum and maximum eigenvalues of the matrix 
M = X^X, and X is the measurement matrix [9], [12], [13]. To mention another example, 
several stability problems in physics, in complex networks and in complex ecosystems are related 
to the probability that all eigenvalues of a random symmetric matrix (for instance with Gaussian 
entries) are negative [14]-[17]. This probability is also important in mathematics, as it is related 
to the expected number of minima in random polynomials [18]. 

Owing to the difficulties in computing the exact marginal distributions of eigenvalues, asymptotic 
formulas for matrices with large dimensions are often used as approximations. One important 
example is the Wigner semicircular law, giving the asymptotic distribution of a randomly picked 
eigenvalue of a symmetric/hermitian random matrix M having zero mean independent, identically 
distributed (i.i.d.) entries above the diagonal. The law applies to a wide range of distributions 
for the entries [19], and in particular when M = X + X^ and the elements of X are zero mean 
(real or complex) i.i.d. Gaussian. In this situation, the symmetric matrix M = X + X^ belongs 
to the Gaussian Orthogonal Ensemble (GOE) and to the Gaussian Unitary Ensemble (GUE), for 
the real and complex cases, respectively. 

Another key result is the Marchenko-Pastur law, giving the asymptotic distribution for one 
randomly picked eigenvalue of the random matrix M = X^X. This is related to the sample 
covariance matrix, and therefore of primary importance in statistics and signal processing. The 
limiting Marchenko-Pastur law applies to a wide class of distributions for the entries of X, 
including the case when X has i.i.d. Gaussian entries, and thus M = X^X is a white Wishart 
matrix [19]. 

The limiting value and the limiting distribution of the extreme eigenvalues (largest or smallest) 
have also been studied intensively [19], [20]. The Tracy-Widom laws give the asymptotic dis¬ 
tribution of the extreme eigenvalues around the limiting values [7], [20]-[27]. Earge deviation 
methods are used in [16], [28], [29] to derive the asymptotic behavior of the distribution of the 
largest eigenvalue far from its mean value. 


Submitted 2015 


DRAFT 



VOL. X, NO. Y, MONTH 2015 


3 


For non-asymptotic analysis, where the random matrices have finite dimensions, deriving the 
distribution of eigenvalues is generally difficult, especially for the real Wishart and GOE. Re¬ 
cently, the exact distribution of the largest eigenvalue, as well as efficient recursive methods for 
its numerical computation, has been found for real white Wishart, multivariate Beta (also known 
as double Wishart), and GOE matrices [30], [31]. 

In this paper we give new expressions and efficient recursive methods for the evaluation of 
the probability ip{a,h) = P{a < Amin(M), Amax(M) < h) that all eigenvalues are within an 
arbitrary interval [a, 6], when M is a real finite dimensional white Wishart, double Wishart, or 
Gaussian symmetric matrix. Eor completeness we provide also the results for complex Wishart 
(with arbitrary covariance), complex double Wishart, and GUE, by specializing [32, Th. 7]. The 
marginal cumulative distribution of the smallest eigenvalue and of the largest eigenvalue can be 
seen as the particular cases 1 — ip{a, cx)) and ip{—oo, h), respectively. 

We then derive simple and accurate approximations to ilj{a,b) based on gamma incomplete 
functions and valid for large matrices, and prove that the probability that all eigenvalues are 
within the limiting spectral support (given by the Marchenko-Pastur and the semicircle laws) 
tends for large dimensions to the universal values 0.6921 and 0.9397 for the real and complex 
cases, respectively. 

Throughout the paper we indicate with r(.) the gamma function, with 7 (a; x, y) = 
the generalized incomplete gamma function, with P{a,x) = j^7(a;0,a;) the regularized lower 
incomplete gamma function, with P(a; x, y) = P~^e~^dt = P{a, y) — P{a, x) the general¬ 

ized regularized incomplete gamma function, with B{a, b) the beta function, with B (x, y, a, b) = 
—ty~^dt the incomplete beta function [33, Ch. 6], with ()1 transposition and complex 
conjugation, and with | ■ | or det(-) the determinant. When possible we use capital letters for 
random variables, and bold for vectors and matrices. We say that a random variable Z has a 
standard complex Gaussian distribution (denoted CJ\f{0, 1)) if Z = {Zi +iZ 2 ), where Zi and Z 2 
are i.i.d. real Gaussian A/’(0,1/2). A complex random vector X is Gaussian circularly symmetric 
if its probability distribution function (p.d.f.) has the form /(x) oc exp (—xlS“^x), where S is 
the covariance matrix. Note that this implies that X is zero mean. When S = I the entries of 
X are i.i.d. CJ\f{0, 1). 
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II. Real Wishart and Gaussian symmetric matrices 

A. Wishart real matrices 

Assume a Gaussian real nxN matrix X with i.i.d. columns, each with zero mean and covariance 
S. The real matrix M = XX^ is called Wishart, and its distribution indicated as >Vn(iV, S). 
When S oc I the matrix is called white or uncorrelated Wishart. 

Denoting Tm{a) = — {i — l)/2), the joint p.d.f. of the (real) ordered 

eigenvalues Ai > A 2 ... > > 0 of the real Wishart matrix M ~ Wn(iV, I) is given by 

[1], [34] 

n n 

/(xi, ...,Xn) = K JJ JJ {xi - Xj) (1) 

2=1 i<j 

where Xi> X 2 > ■ ■ ■ > Xn>^, a = {N — n — l)/2, and K is a normalizing constant given by 

^ ^ 2^^/^Tn{N/2)Tn{n/2) ' 

B. Gaussian symmetric real (GOE) matrices 

The Gaussian Orthogonal Ensemble is constituted by the real n x n symmetric matrices whose 
entries are i.i.d. Gaussian Af(0,1/2) on the upper-triangle, and i.i.d. A/'(0,1) on the diagonal 
[20]. The joint p.d.f. of the eigenvalues for GOE is [3], [20] 

n n 

f(xi, ...,Xn) = Kgoe JJ JJ (Xi - Xj) (3) 

2=1 i<j 

where Ti > 0:2 > ■ ■ ■ > and Kgoe = YYe=i r(^/2)]“^ is a normalizing constant. Note 
that here the eigenvalues are distributed over all the reals. 

C. Real multivariate beta (double Wishart real) matrices 

Eet X, Y denote two independent real Gaussian p x m and p x n matrices with m, n > p, 
each constituted by zero mean i.i.d. columns with common covariance. Multivariate analysis of 
variance (MANOVA) is based on the statistic of the eigenvalues of (A + B)^^B (beta matrix), 
where A = XX^ and B = YY^ are independent Wishart matrices. These eigenvalues are clearly 
related to the eigenvalues of A^^B (double Wishart or multivariate beta). 
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The joint distribution of s non-null eigenvalues of a multivariate real beta matrix in the null case 
can be written in the form [2, page 112], [1, page 331], 

S S 

f{xi,...,Xs) =Kmb Wx'^{l - XiY^{xi - Xj) (4) 

2=1 i<j 

where 1 > xi > 0:2 ■ ■ ■ > > 0, and Kmb is a normalizing constant given by 

Y' ^ i+2m~\-2n-\-s-\-2 ^ 


Kmb — 


(5) 


iir(|)r(^)r(^) 

With the notation introduced above, this is the distribution of the eigenvalues of (A + B)^^B 
with parameters s = p, m = (n — p — l)/2, n = (m — p — l)/2 . The marginal distribution of 
the largest eigenvalue is of basic importance in testing hypotheses and constructing confidence 
regions in MANOVA according to the Roy’s largest root criterion [1, page 333], [31]. 

D. The function h) for real Wishart matrices 

The following is a new theorem for real white Wishart matrices. 

Theorem 1. The probability that all non-zero eigenvalues of the real Wishart matrix M ~ 
>V„(A, I) are within the interval [a, 6] C [0, cx)) is 

f{a, b) = Pr(a < A^i„(M), A^ax(M) <h) = K' ^\A{a,b)\ (6) 

with the constant 

n 

K' = K r (a + f) . 

t=i 

In (6), when n is even the elements of the n x n skew-symmetric matrix A(a, b) are^ 

2 


O-iJ — 


K [oij,-] + P (aj, - 


a b 

P ( 2 ’ 2 


r(aj 


X 


a+i—l—x 


'a/2 


e ^P{aj,x)dx (7) 


for i,j = l,...,n, where ai = a ^ = (A — n — l)/2 + T 

When n is odd, the elements of the (n + 1) x (n + 1) skew-symmetric matrix A(a, b) are as in 


(7), with the additional elements 

/ a b\ 



*^2,n+l 

^("'’2’ 2 ) 

i = 1,... ,n 


^n+lj 

^j,n+l 

j = l,...,n 

(8) 

<^r2+l,n+l 

= 0. 




'Note that skew symmetry implies Uij = —aj,i and ai^i = 0. 
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Moreover, the elements a*j can be computed iteratively, without numerical integration or series 
expansion, starting from 0 *,* = 0 with the iteration 




r(aj + aj) 

T{aj + l)r(ai)2"'+“i-i 


P{ai + aj] a, b) 


g{aj,a/2) + g{aj,h/2) / a 

r(a, + l) P’2’2y' 

(9) 


for j = i,... ,n — I, with g{a, x) = a;“e 


Proof: First we observe that an identity related to Vandermonde matriees gives det [{yl } = 


- Vi)- 

Then, for arbitrary eonstants 7 ^ 0 we write 


f{xn,...,xi)d:s: = K' 




/ P n n 

■ n 9 '* Tl(yj - yi)^y 


“/2 <yi<"-<yn< b/2 


det [{$*(%•)}] dy. 


( 10 ) 


^ ] ■ J 

«/2 <yi<---<yn< b/2 

with ^i{y) = 7 i y^+^-^e~y and K' = 11”=! 77^- evaluate this integral we reeall 

that for a generie m x m matrix #(w) with elements where the ^fx), i = 1,... ,m 

are generic functions, the following identity holds [35] 

J-J |$(w)|dw = Pf(A) ( 11 ) 

a<wi<...<Wm^b 

where Pf(A) is the Pfaffian, (Pf(A))^ = det A, and the skew-symmetric matrix A is m x m 
for m even, and (m + 1) x (m + 1) for m odd, with 


b pb 


Q-iJ — 



( 12 ) 


sgn(|/ — x)^i{x)^j{y)dxdy i, j = 1,... ,m. 

For m odd the additional elements are ai^m+i = —(im+i,i = Ja^i{x)dx, i = and 

®m+l,m+l 0. 

We then note that, by writing Fi{y) = ^e{x)dx a primitive of ^eix), we can rewrite (12) as 


b r rb 


0-i,j 


<^i{x)(^j{y)dy - / ^i{x)^j{y)dy 


<hi(a;) [Fj{b) + Fj{a) — 2Fj{x)] dx 


dx 


= [F,{b) + F,(a)] [Ffb) - Fi{a)] - 2 / ^fx)F,{x)dx . 


(13) 
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Then, using (11) and (13) in (10) with $i(a;) = * and 7 ^ = l/r(a + i), after some 

manipulations we get ( 6 ) and (7). 


Theorem 1 does not require numerical integration or infinite series. In fact, first we observe that 
for an integer n we have [33, Ch. 6 ] 


P(a + n,x) = P(a, x) 


n—1 


k=0 


^a+fc 

r(a + /c + 1 ) 


(14) 


Therefore, P(a, x) and P(a; x, y) can be written in closed form when a is integer or half-integer, 
starting from P(0,a;) = 1 and P(l/2,a;) = erf^^. Moreover, using the relation P(a + l,a;) = 
P{a,x) — e~^x°‘/T{a + 1) in (7) gives, after simple manipulations, the iteration (9). 


In summary, the probability that all eigenvalues are within the interval [a, h] is simply obtained, 
without any numerical integral, by Algorithm 1. For instance, implementing directly the algorithm 
in Mathematica on a personal computer we obtain the exact value ^(a, 6 ) for n = iV = 500 in 
few seconds. 


E. The function ■^(a, h) for real symmetric Gaussian matrices (GOE) 
The following is a new theorem for GOE matrices. 


Theorem 2. The probability that all eigenvalues of the real GOE matrix M are within the 
interval [a, 6] C (—cx), cx)) is 

f{a, b) = Pr(a < Ami„(M), Amax(M) <h) = \/|A(a,6)| 

with the constant 

n 

K'goe = Kgoe n r (^/ 2 ) • 


£=1 


In (15), when n is even the elements of the nx n skew-symmetric matrix A{a,b) are 


for i,j = 1 ,... ,n. 


P — + P, — 


FI — P 


V2j\ ^\V2^V2j rmj a/72 


-6/72 


p ^"Fj{x)dx (16) 


(15) 
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Algorithm 1 (a, 6) for real Wishart matrices 

Input: n, TV, a, b 

Output: ij{a, b) = Pr {a < Amin(M), Amax(M) < b} 

A = 0 


ai = (TV — n — l)/2 + f 


x) = 

for i = l—)-n — Ido 
for j = T —)■ n — 1 do 




r(a, + aj) 

r(Q;j + l)r(Q;j) 


3 

-P{ai+aj-, a, b) 


g{aj,a/2) + g{aj,b/2) 
T{aj + 1 ) 


end for 
end for 

if n is odd then 

append to A one column according to (8) and a zero row 

end if 


P 


) 2'2 


A = A - A^ 
return K' ^/\A\ 


When n is odd, the elements of the (n + 1) x (n + 1) skew-symmetric matrix A(a, b) are as in 
(16), with the additional elements 



^ f a b \ 





T = 1,..., n 


^n+lj 

^7,72+1 

j = 

(17) 

*^n+l,n+l 

= 0. 




In the previous equation 


F, (y) = 


r(j/2) 


I'y . , 1 

^ rinr* — 


e ^ dx =-sgn^{y)P ( -,y 


3 .2 


(18) 


and Fj{x,y) = Fj{y) - Fj{x). 

Moreover, the elements aij can be computed iteratively, without numerical integration or series 
expansion, starting from 


02,1 = ^ I [erfc{b) - eifc{a)] + (^e + e erfc - ^^rfc 


b 


(19) 


Submitted 2015 


DRAFT 




VOL. X, NO. Y, MONTH 2015 


9 


and using the antisymmetry aj^i = — ajj, together with the iteration 








r(i/2)r(j/2 +1) ‘+^ 

where q{j,x) = x^e~^ . 


{a,b) 


q{j,a/V2) + q{j,b/V2) 
2T{j/2 + 1) 


F- — — 

\V2^V2 


( 20 ) 


Proof: For arbitrary constants 7 * 7 ^ 0, we have 



fx{xn,.. .,xi)dx: 


a <xi<-"<Xn^ b 


K' 

^GOE 


Wlie ^^^W{yj-yi)dy 


^l\^<yi<-"<yn< Vv^ 


i=l 


i<j 


K' 

^GOE 


det[{<l>i(|/j)}]ciy. (21) 


ilV2<yi<---<yn< h/\/2 


with ^fy) = yiy^ and nr=i 7* Then, using (11) in (21) with 

<hj(a;) = 7 ia;*“^e“^^ and 7 * = l/r(j/ 2 ), after some manipulations we get (15). 

Theorem 2 is amenable to easy evaluation, without numerical integration or infinite series. In 
fact, first we observe that, substituting (14) in (18) we have 


Fj+2{y) 


Fj{y)-y^e F 


1 

2 r(j 72 + i)- 


Using also the relation 

^ 2-^-F^sgn{f3)^T Q) P (|,27^) 

in (16) gives, after some manipulations, the iteration ( 20 ). ■ 

To build iteratively the upper half of the skew-symmetric matrix A{a,b) we just need (20) 
and the first two diagonals a* j and Oj j+i. The first diagonal is clearly identically zero due to 
skew-symmetry, giving 0^7 = 0. The odd first diagonal Oi^i+i can be obtained by a zig-zag 
iteration 

(a) (b) (a) {b) (a) 

fll,2 -^ 02,1 -^ 02,3 -^ 03^2 -^ 03^4 -?■ 04^3 ■ ■ ■ 


where steps (a) use skew-symmetry, and steps ( 6 ) use (20). The element 01^2 is directly obtained 
in closed form from (16). 

In summary, the probability that all eigenvalues are within the interval [a, b] is simply obtained, 
without any numerical integral, by Algorithm 2. 

The algorithm can be used to evaluate numerically or symbolically f{a,b). Evaluating numer¬ 
ically ilj{a,b) for an arbitrary interval [a,b] requires few seconds for matrices of dimensions 
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Algorithm 2 h) for real Gaussian symmetric matrices (GOE) 

Input: n, TV, a, b 

Output: ^|J{a, b) = Pr {a < Amin(M), Amax(M) < b} 

A = 0 

ai ,2 = — \ |v^ [erfc(6) — erfc(a)] + 
for i = l— )-n — 2do 
for j = i —)■ n — 2 do 

derive 0*^+2 from aij using (20) 

end for 

derive aj+i,j +2 from using (20) 

O'i+lji 0 

end for 

if n is odd then 

append to A one column according to (17) and a zero row 

end if 

A = A - A^ 

return K'^oe VW\ 



n = 500. The exact expression of ^(a, b) can be also derived symbolically in closed form. Some 
examples for the probability that all eigenvalues are positive, obtained from Algorithm 2, are: 


n = 1 
n = 2 

n = 3 


^(-00, 0 ) 
^(-00, 0 ) 
^(-00, 0 ) 


1 

71 - 2^2 

Att 
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n = 4 
n = 5 

n = 10 


^(-oo,0) 

^(-oo,0) 


y 1(9 - 4v^) (-16 - 4v^ + Ttt) 
SGtt 

—8 — + Stt 

24^r 


^(-oo,0) 


(44217 - 27392^2) 

1833775104007r2 

(^432799744 + 6251520^2 - (^278413220 + 1989925^2) tt + 447699007r2 


The expressions for n = 1, 2 and 3 were already known as reported in [16]. 

In [16] the following asymptotic bound is also derived 

cx),0) ^ (-22) 

Higher order corrections for large n have been provided in [29, eq. (19)], not reported here for 
space reason. By comparing the exact value of ^(—cx),0) with the approximation above, we 
found that the error is exponential in n, and well approximated as a factor 10“”/®. We found 
therefore that an improved approximation for ^(—cx), 0) is 

^7(-00, 0) ^ e“”" l°(3)/4-nln(10)/6 _ (23) 

Some values are reported in Table I. 


TABLE I 

Probability ^{—oo, 0) that all eigenvalues are negative, GOE. 


n 

exact (Alg. 2) 

approx. (22) 

approx. [29, eq. (19)] 

approx. (23) 

2 

0.146 

0.333 

0.322 

0.155 

5 

1.40E-4 

1.04E-3 

1.91E-3 

1.53E-4 

10 

2.27E-14 

1.18E-12 

1.23E-12 

2.54E-14 

50 

2.43E-307 

6.30E-299 

3.31E-304 

2.92E-307 

100 

2.72E-1210 

1.57E-1193 

1.49E-1206 

3.39E-1210 

500 

2.85E-29904 

8.35E-29821 

3.88E-29899 

3.87E-29904 


Similar considerations can be done for Wishart matrices, for which approximations for large 
deviation behavior of the largest eigenvalue are available [36]. In particular, [36, eq. (4)] is an 
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approximation for ■^(0, n) for Wishart matrices M Wn(n,I). In Table II we eompare the exaet 
results and the approximation for some values of n. 

TABLE II 

Probability i/)( 0, n), real Wishart, N = n . 


n 

exact (Alg. 1) 

approx. [36, eq. (4)] 

2 

0.315 

0.491 

5 

3.71E-3 

L18E-2 

10 

L90E-9 

L95E-8 

50 

L70E-198 

L81E-193 

100 

10.2E-781 

L07E-771 

500 

7.33E-19325 

6.22E-19275 


F. The function for real multivariate Beta (double Wishart) matrices 

The following is a new Theorem for multivariate real beta matriees in the null ease. 


Theorem 3. The probability that all eigenvalues of a real multivariate beta matrix M are within 
the interval [a, b] C [0,1] is 


fia, b) = Pr(a < A^in(M), A^ax(M) <b) = K' s/\A{a,b)\ 


(24) 


with the constant 


^MB — AImB 


r (m + 


(25) 


r (m + f + n + 1 ) ’ 

In (24), when s is even the elements of the s x s skew-symmetric matrix A(a,b) are 

cii,j = kikj [13(0, a; m + j, n + 1) + 13(0, 6; m + j, n + 1)] 13(a, 6; m + i, n + 1) + 

— 2kikj f — x)^B{0,x;m + j,n + l)dx (26) 

J a 

for 1 , j = 1, ..., s, where ki = T(m + n + f + l)/r(m + i). 
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When s is odd, the elements of the (s + 1) x (s + 1) skew-symmetric matrix A(a, b) are as in 
(26), with the additional elements 


0'i,s+l 

= kiB{a,b]mi,n1) 






7 = 1,- 

..,s 

(27) 


= 0 





Moreover, the elements Oij can be computed iteratively, without numerical integration or infinite 
series expansion, starting from Oi^i = 0 with the iteration 


ij+i = ai^j-ki[gj+i{a) + gj+i{h)]B{a,h]m+i,n+l) + 


m + n + j + 1 

for j = i,... ,s — 1, with ge{x) = {m + n + i). 


‘^kikj+i B{a,b]2m-\-i-\-j,2n-\-2) 


(28) 


Proof: Similarly to the previous, the proof leading to (26) uses (11), (12) and (13), with 
4>i{x) = — xY and F^y) = Jq = keB{0, x; m + f, n + 1). 

For the iterative derivation of the elements aij, we use the property [31] 

B{0, x; a + 1, &) = —^ig(0, x; a,b) . (29) 

a + 0 a + 0 

which produces 

Fj+i (x) = Fj (x) - gj+i (x). (30) 


In summary, the probability that all eigenvalues are within the interval [a, b] is given by Algo¬ 
rithm 3. Implementing directly the algorithm in Mathematica on a personal computer, we obtain 

for example the exact distribution of the largest eigenvalue in less than 0.1 seconds for all tables 
in [37], [1, Table B.4] and [25, Table 1], 
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Algorithm 3 (a, 6) for real multivariate beta matriees 

Input: s, m, n, a, b 

Output: ij{a, b) = Pr {a < Amin(M), Amax(M) < b} 

A = 0 


ge{x) = ^{1 — x)^~^^k£/{m + n + i) 

kc = r{m + n + i + l)/r(m + i) 


for i = l—)-s — Ido 
for j = i — )■ s — 1 do 

a*j+i = ttij - ki [gj+i{a) + gj+i{b)] B{a, b-,m + i,n + 1) + 
i + j, 2n + 2) 

end for 
end for 

if s is odd then 

append to A one eolumn aeeording to (27) and a zero row 

end if 


2kikj-^\ 
m + n + j + 


-B{a, b] 2m + 


A = A - A^ 

return K'^^^/\A\ 


III. Complex uncorrelated and correlated Wishart and hermitian Gaussian 

MATRICES 

The analysis for eomplex random matriees is easier than for the real ease, and in faet some 
important results are known sinee many years for eomplex multivariate Beta matriees and 
for uneorrelated eomplex Wishart [38]. A general methodology whieh ean be applied also to 
eorrelated eomplex Wishart (i.e., with eovarianee matrix not proportional to the identity matrix) 
is given in [32] and here speeialized to provide ip{a,b) in several situations. 

A. Complex Wishart 

Assume a Gaussian eomplex n x N matrix X with i.i.d. eolumns, eaeh eireularly symmetrie 
with eovarianee S = I, and N > n. Denoting Tn{m) = ~ 0-’ the joint 

p.d.f. of the (real) ordered eigenvalues Ai > A 2 ... > An > 0 of the eomplex Wishart matrix 
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M = XX^ ~ CWn(-/V, I) (identity covariance) is well known to be [5], [7], [8] 

n n 

f{xi, ...,Xn) = K JJ JJ (Xi - Xjf (31) 

2=1 i<j 

where Xi > X 2 > ■ ■ ■ > Xn > 0 and 7^ is a normalizing constant given by 

n 

l/K = IK^ -i)\{n-i)\ (32) 

2=1 

Assume now a Gaussian complex nx N matrix X with i.i.d. columns, each circularly symmetric 
with covariance S, and N > n. The joint distribution of the ordered eigenvalues of M = XX^ ~ 
CWn{N, S) has firstly been found in [8] as follows. 

Lemma 1. Let M ~ C’Wn(iV, S) be a complex Wishart matrix, N > n. Denote ai > a 2 > 

... > > 0 the ordered eigenvalues of S. Then, the joint p.d.f of the ordered eigenvalues of 

M is 

n n 

f{xi,... ,Xn) = K-e\E (x, cr)| ■ ]^(a;i - xf ■ (33) 

i<j j=l 

where E (x, cr) = and 

n n 

1/A'e = Iff'"' - n - 0!. (34) 

i<j 2=1 

Proof: See [8]. ■ 

The analysis in [8] has been extended to the case where S has eigenvalues of arbitrary multiplicity 
and to the marginal eigenvalues distribution in [32], [39], [40]. 

In particular, when S is spiked with ui > (T 2 = <^3 = <^4 = ■ ■ ■ = we have the following 
result. 


Lemma 2. Let M ~ C’W„(A^, S) be a complex Wishart matrix, N > n. Denote ai > a 2 = 

... = (Tn > 0 the ordered eigenvalues of S (spiked covariance matrix). Then, the joint p.d.f of 
the ordered eigenvalues of M is 


f(xi, ...,Xn) = Ki |E(x,cr)| ■ JJ(a;i - xf ■ JJaif 

i<j j=l 


where E (x, cr) has elements 


e 








7 = 1 

j = 2 ,..., n 


(35) 


(36) 
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l/Ki = - a2T~^ - *)! n 

i=l 1=2 

Proof: This is a particular case of [39, Lemma 6]. ■ 

Below we report ■?/’(«, h) for complex Wishart matrices. 

Theorem 4. For complex Wishart matrices M ~ CWn(iV, S), N > n, the probability that all 
eigenvalues are within [a, h] C [0, cx)) is given below, depending on the covariance S. 

1) For the uncorrelated complex Wishart matrix M ~ CWn{N, I).- 

Pr {a < Amm(M), Amax(M) <b} = K |A(a, b)\ (38) 

where the elements of the n x n matrix A (a, b) are 

Qij = / = 7 (A + n — i — j + 1; a, 6) (39) 

J a 

and K is given in (32). 

2) For the correlated complex Wishart matrix M ~ CWn(A, S) where S has distinct 
eigenvalues ai > a 2 > ■ ■ ■ > an: 

Pr{a < Amin(M), Amax(M) < 6} = As |A(a,6)| (40) 

where the elements of the n x n matrix A (a, b) are 

aij = ^ + (41) 

and As is given in (34). 

3) For the correlated complex Wishart matrix M ~ CWn{N, S) with a spiked covariance S 
having eigenvalues ai > a 2 = a^ = a^ = ■ ■ ■ = an: 

Pr{a < Amin(M), Amax(M) < 6} = Ai |A(a,6)| (42) 


where the elements of the n x n matrix A (a, b) are 




fN-^e-^/^^dt = - * + 1 ; —, — 1 

V 0-1 0-1 / 


J = 1 


fN+n-i-je-F-^dt = af (N + n-i-j + 1;—,—'] j = 2,..., n 

I V (^2 (^ 2 / 
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and Ki is given in (37). 

Proof: This Theorem can be obtained by specializing [32, Theorem 7]. More precisely, 
we first rewrite the p.d.f.’s in (31), (33), and (35) as product of determinants of two matrices, 
by expressing Y[i<j{^i ~ ^j) the determinant of a Vandermonde matrix. Then, applying [32, 
Theorem 7], after some simplifications we get the Theorem. ■ 

Note that the uncorrelated case 1) in the previous theorem can be seen as an extension of [38, 
eq. (6)]. 

We remark that approximations and asymptotics for spiked Wishart have also been studied in 
recent literature (see e.g. [7], [19], [41]). 

B. Complex multivariate beta (double Wishart complex) matrices 

When X, Y are two independent complex Gaussian, the analogous of (4) is the complex multi¬ 
variate beta, where the joint distribution of the eigenvalues is [38] 

S S 

f{xi,...,Xs) = Kmb Wxf{l - Xi)"^ -^{Xi - Xjf (44) 

i=l i<j 

with 1 > xi > a ;2 ■ ■ ■ > > 0, and 

_TT r(m + n + s + i) 

"Li r (*) T (i + rn)T (i + n) ' 

Therefore, by applying [8, Corollary 2] we have for a complex multivariate Beta matrix M 

f{a, b) = Pr {a < Amin(M), Amax(M) <h} = Kmb | A(a, h) \ (45) 

where the elements of the s x s matrix A (a, h) are: 

ttij = B(a, b] m + i + j — 1, n + 1) (46) 

for i, j = 1,..., s. 

Note that (45) can be seen as an extension of [38, eq. (3)]. 

C. Hermitian random Gaussian matrices (GUE) 

The Gaussian Unitary Ensemble (GUE) is composed of complex Hermitian random matrices 
with i.i.d. CA/’(0,1/2) entries on the upper-triangle, and A/'(0,1/2) on the diagonal [20]. 

The following Theorem applies to the GUE. 
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Theorem 5. The probability that all eigenvalues of a n x n GUE matrix M are within the 
interval [a, &] C (—cx), cx)) is 

f{a,b) = Pr{a < Amin(M), Amax(M) < b} = |A(a, 6)| (47) 


where the elements of the n x n matrix A (a, b) are 

pb 

aij= / dt 


P ( ^^^ 4 —bA sgn{by^^ ^ -P ( ^^- 4 —■’ 


(48) 

i+jr-ll 


I f t+J-l 

2V 2 yrv 2 ’“y ' V 2 

and Kgue = \Xi=i ^ normalizing constant. 

Proof: As for the complex white Wishart, this theorem for GUE is easily derived from 
known results. In fact, the joint distribution of the ordered eigenvalues can be written as [20] 

n n 

f{xi, ...,Xn) = Kgue JJ {xi - Xjf JJ . (49) 

i<j 2=1 

Then, by using for example [8, Corollary 2] with '^fxf = ^fxj) = x^~^,^{x) = we get 
immediately the result. ■ 


IV. Asymptotics and approximations 

In this section we study f{a, b) for large white Wishart and Gaussian matrices. To this aim we 
make the following three observations. 

1) The statistical dependence between the largest and the smallest eigenvalues passes through 
the intermediate n — 2 eigenvalues. Consequently, in the limit for n —)■ cx the largest 
eigenvalue and the smallest eigenvalue are independent. Thus, for large matrix sizes we 
have 

Pr (a < Amin(M), Amax(M) < 6) Ri Pr (a < Amin(M)) Pr (Amax(M) < b) (50) 

which is like to say ip{a, b) ip{a, oo)f{—oo, b). 

2) The distribution of the smallest and largest eigenvalues of white Wishart and Gaussian 
matrices tend to a properly scaled and shifted Tracy-Widom distribution [7], [20]-[23], 
[26], [27], [42]. 
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More precisely, focusing for example on white Wishart matrices M ~ >Vp(m,I) or M 
CWp(m,I), when m,p — )■ cx) and m/p — )■ 7 G [0, ex] 

-^max (M) Pmp V 


a, 




(51) 


mp 


where denotes the Tracy-Widom random variable of order /) whose cumulative 

distribution function (CDF) will be indicated with F^^x), and^ 


Pmp — T y/p) (^mp — y/JFnp [l/ y/P T 1 


1/3 


(52) 


In the previous expressions /? = 1 and /3 = 2 for real and complex matrices, respectively. 
Similarly, for the smallest eigenvalue we can use the results in [26] or [27]. For example, 
in [26] it is shown that when —)■ x and m/p —)■ 7 G (1, x), 

-^min(M) Prrip V ^ 


a. 


rwp 


mp 


with 


Pmp = [Vm 


-VpY [ 1 /v^- 1 /v^] 


1/3 


(53) 


(54) 


Note that the variance of the smallest eigenvalue is smaller than that of the largest 
eigenvalue. 

For the Gaussian ensemble we have^ for n —)■ x [20]-[22] 

-^max (M) — fin V 


(Jri 


'^min(h/[) Pn ® 


with 


fin = 2ao^n 




Pn ~ Pri 




TWf 


(Jn = dn = Cio{n) 


(55) 

(56) 

(57) 


3) The Tracy-Widom distribution can be accurately approximated by a scaled and shifted 
gamma distribution 

rW /3 ^ T{k, e)-a (58) 


^For small sizes a better approximation is obtained by using slightly different values, like m — 1/2 and p — 1/2 for the real 
case, instead of m,p. However, since ip{a,b) can be computed exactly as shown in Section II and III, asymptotic expressions 
are interesting only for large dimensions for which these corrections are irrelevant. 

^For GOE and GUE the smallest and largest eigenvalues have symmetrical distributions. 
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TABLE III 

Parameters eor approximating TW/3 with r[fc, 6] — a. 



TWi 

TTVa 

TWi 

k 

46.446 

79.6595 

146.021 

e 

0.186054 

0.101037 

0.0595445 

a 

9.84801 

9.81961 

11.0016 


where a is a constant, and T{k,9) denotes a gamma random variable (r.v.) with shape 
parameter k and scale parameter 9 [30]. Thus the CDF of TW /3 is accurately approximated 
by an incomplete gamma function as: 

Pr(rW /3 <x) = F/ 3 {x) ~ P ^ 

where (a;)+ = max{0, x} denotes the positive part. The parameters k, 9, a are reported in 
Table III [30]. 



Thus, putting the gamma approximation in (51), (53) we have for white Wishart, GOE and GUE 
matrices 


Pr (Amax(M) < 6) —)■ P/3 


Pr (Amin(M) > a) —)■ P/3 


b — fj, 


a 

a — jjL^ 

CT' 


~ P A;, 


(« + (&- /i)/a)' 


~ P ( A;, 


9 


{a — {a — fj, )/a )’ 
9 


(60) 


(61) 


where , a are given by (52), (54), (57), and the parameters k, 9, a are given in Table III. 

These can be used in (50) to give 


b) ^ P i k, 


{a + {b- 
9 


-\ + 


P A;, 


{a — {a — jj, )/a ) 

T~ 


(62) 


Einally, we observe that (62) can be used not only for white Wishart and Gaussian symmet- 
ric/hermitian matrices, but for a much wider class of matrices, due to the universality of the 
Tracy-Widom laws for the smallest and largest eigenvalues of large random matrices [24], [26], 
[43]. 


V. Probability the all eigenvalues are within the support of the limiting 
Marchenko-Pastur (MP) and Wigner spectral distribution 

Under quite general conditions, for a matrix M = XX^ where X is (p x m) with i.i.d. entries 
with zero mean and variance = 1, the Marchenko-Pastur law gives the asymptotic p.d.f. of 
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an unordered"^ eigenvalue A = A(M) for large p,m with a fixed ratio p/m < 1 as^ 


f{X) = { 27ip\ 
0 


{b — A)(A — a) a < X <b 
otherwise 


(63) 


where a = [y/m — a/p)^ and b = [y/m + a/p)^, and [a, b] is the support of the MP law [19], 
[44], 


Also, for increasing p, m it has been proved that the the minimum and maximum eigenvalues 
converge to the edges of the MP law Amin(M) —)■ a and Amax(M) —)■ b [19]. Note that when the 
entries of X are Gaussian the matrix M is white Wishart. 

Similarly, under quite general conditions, for a Wigner matrix M = X + X^ where X is (p x p) 
with i.i.d. entries with zero mean and variance = 1/4, the Wigner semicircle law gives the 
asymptotic p.d.f. of an unordered eigenvalue A = A(M) for large p as 


/(A) 


— V2p-A2 |A| < 

Tip 

0 otherwise 


(64) 


where [—y/2p^ is the support of the semicircle law. Again, for increasing p it has been 
proved that the the minimum and maximum eigenvalues converge to the edges of the semicircle 
support Amin(M) -)■ -^/2p and Amax(M) -)■ ^2/) [19]. 


So, we could be tempted to think that for increasing matrix sizes all eigenvalues are within the 
MP or semicircle supports with probability tending to one. However, this is not the case, as 
proved in the following Theorem. 


Theorem 6. For increasingly large matrices, the probability that all eigenvalues of Wishart and 
Gaussian hermitian matrices are within the MP and semicircle supports tends to F/(0) ~ 0.6921 
and Tf(0) ~ 0.9397/or the real and complex cases, respectively. 

More precisely, we have the following results. 

1) Let M ~ Wp{m, I) be a real Wishart matrix with m > p. When m,p ^ oo and m/p —)■ 
7 G (1, oo), the probability that all eigenvalues are within the Marchenko-Pastur support 


''This can be seen as the distribution of a randomly picked eigenvalue, or of the arithmetic mean of all eigenvalues. In physics 
literature, this is related to the fraction of eigenvalues below a given value (spectral distribution). 

^If p/m > 1 the matrix has p — m zero eigenvalues, so the distribution has an additional point of mass 1 — m/p in 0. 
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is 

%l) {{^/m - a/p)^, (a/th + a/p)^) -^1 (0) ~ 0.6921. (65) 

2) Let M ~ C>Vp(m, I) be a complex Wishart matrix with m > p. When m,p —)■ cx) and 
m/p —)■ 7 G (1, oo), the probability that all eigenvalues are within the Marchenko-Pastur 
support is 

%l) ((aA« - a/p)^, + a/p)^) -^2^(0) ^ 0.9397. (66) 

3) Let yi be a {p x p) real symmetric GOE matrix. When p ^ oo the probability that all 
eigenvalues are within the semicircle support is 

i) ^ Fi^(O) ~ 0.6921. (67) 

4) Let yi be a {p X p) complex symmetric GUE matrix. When p ^ oo the probability that 
all eigenvalues are within the semicircle support is 

ilj ^ F|(0) ~ 0.9397. (68) 

Proof: Let us consider the Wishart case. The statistical dependence between the largest 
and the smallest eigenvalues passes through the intermediate p — 2 eigenvalues, and in the limit 
for p —)■ CX) the largest and smallest eigenvalues are thus independent. Then, using (50) and 
observing that the MP edges are the constants pmp and appearing in (51) and (53), we get 
the results. For GOE/GUE the limiting distribution of the extreme eigenvalues is still a shifted 
(of an amount equal to the circular law edges) and scaled TWp, so the same reasoning leads 
to (67) and (68). ■ 

This theorem is valid not only for matrices derived from Gaussian measurements, but for the 
much wider class of matrices for which (51) (53) applies [24], [26], [43]. 

Note that the gamma approximation (59) for the Tracy-Widom law gives (0) ^ (/c, |) = 

0.83122 _ Q gg^ ^2(0) _ p2 ^ 0.969452 ~ 0.9398. 

Exact values for ^(a, h) for finite dimension real Wishart matrices, as obtained by Algorithm 1, 
are reported in Table IV, together with the asymptotic value. 

Einally, since for Wishart matrices the smallest and largest eigenvalues have different variances, 
asymptotically we can leave the same probability at the left and right side by moving t a” ^ on 
the left and t ams on the right of the limiting MP support. In fact, we have 
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Pr (Amax(M) > Hms + t (Jms) -)■ 1 - Fp(t) (69) 

Pr (Amin(M) < - t a“,) ^ 1 - Fp{t) (70) 

and thus 

^ '^msi + t <7ms) ^ (1 Ffjit)) (71) 

with, e.g., Fi(0) ~ 83%, Fi(l) ~ 95%, Fi(2) ~ 99%, Fi(3) ~ 99.8%. 


TABLE IV 

Probability V'(a, b) that all eigenvalues of a real Wishart matrix are within the Marchenko-Pastur 
EDGES FOR p/m = 2/3,1/2,1/5,1/10. Numerical values for finite p obtained by Algorithm 1, and EORp = oo 

BY Theorem 6. 


p 


b) 



p/m = 2/3 

p/m = 1/2 

p/m = 1/5 

p/m = 1/10 

10 

0.7678 

0.7645 

0.7625 

0.7624 

20 

0.7499 

0.7483 

0.7476 

0.7477 

50 

0.7332 

0.7326 

0.7327 

0.7329 

100 

0.7239 

0.7238 

0.7242 

0.7244 

200 

0.7169 

0.7171 

0.7175 

0.7177 

500 

0.7101 

0.7103 

0.7108 

0.7109 

oo 

0.6921 

0.6921 

0.6921 

0.6921 


VI. Application to compressed sensing 
Assume that we want to solve the system 

y = Ax (72) 

where y G MF and A G are known, the number of equations is m < n, and x G M” is the 
unknown. Since m < n we can think of y as a compressed version of x. Without other constraints 
the system is underdetermined, and there are infinitely many distinct solutions x satisfying (72). 
If we assume that at most s elements of x are non-zero (i.e., the vector is s-sparse), and s < m, 
than there is only one solution (the right one) to (72), provided that all possible submatrices 
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consisting of 2s columns of A are maximum rank (2s). However, even when this condition is 
satisfied, finding the solution of (72) subject to | |x| |o < s, where the fo-norm 11 ■ | |o is the number 
of non-zero elements, is computationally prohibitive. A computationally much easier problem 
is to find a fi-norm minimization solution x = {argminx ||x||i : y = Ax}. In [9] it is proved 
that, under some more strict conditions on A, the solution provided by fi-norm minimization is 
the same as that of the fo-norm minimization. More precisely, for integer s define the isometry 
constant of a matrix A as the smallest number 5s = ()s(A) such that [9] 

(1 - ()s)||x ||2 < ||Ax ||2 < (1 + 5s)||x||2 (73) 

holds for all s-sparse vectors x. The possibility to use ii minimization instead of the impractical 
fo minimization is, for a given matrix A, related to the restricted isometry constant [9]. For 
example, in [45] it is shown that the Iq and the ii solutions are coincident for s-sparse vectors 
X if 5, < 0.307. 

Then, the next question is how to design a matrix A with a prescribed isometry constant. One 
possible way to design A consists simply in randomly generating its entries according to some 
statistical distribution. In this case, for given m, n, s and 5s the inequalities (73) are satisfied 
with some probability. The target here is to find a way to generate A such that, for example, for 
given m,n,s, the probability Pr((5s(A) < 0.307) is close to one. When the measurement matrix 
A has entries randomly generated according to a N{0,l/m) distribution, this probability can 
be bounded starting from the probability Pr(a < < b), where A^ is 

a m X s Gaussian random matrix with N{0,l/m) i.i.d. entries [9, Sec. III]. In [9], deviation 
bounds for the largest and smallest eigenvalues of A^A^ are obtained, using the concentration 
inequality, as 

Pr (^y/A„ax(AJA,) > 1 + yZ + 0 ( 1 ) + 

Pr (^y/A^in(A^A,) < 1 - + 0(1) - 

where t > 0 and o(l) is a small term tending to zero as m increases. In our notation, and 
neglecting o(l), the previous bounds can be rewritten 

Pr (A^ax(M) > +V^ + ty/^) ( 74 ^ 

Pr (^Amin(M) < (a/th - y/s - ty/mY^ < (75) 
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where M = mAjA^ ~ Ws(m,I) and Amax(M) = mAmax(AfA 5 ). In Fig. 1 and Fig. 2 these 
bounds are compared with the exact results given by Algorithm 1 and with the simple gamma 
approximations (60), (61), for some values on s,m. It can be noted that the concentration 
inequality bounds (74), (75) are quite loose. 

VII. Conclusions 

Iterative algorithms have been found to evaluate in few seconds the exact value of the proba¬ 
bility that all eigenvalues lie within an arbitrary interval [a,b], for quite large (e.g. 500 x 500) 
real white Wishart, complex Wishart with arbitrary correlation, double Wishart, and Gaussian 
symmetric/hermitian matrices. 

For large matrices, simple approximations based on shifted incomplete gamma functions have 
been proposed, and it is proved that the probability that all eigenvalues are within the limiting 
support is 0.6921 for real white Wishart and GOE, and 0.9397 for complex white Wishart and 
GUE. 

Eor example, we analyzed the probability that all eigenvalues are negative for GOE, of interest 
in complex ecosystems and physics, and compared the new expressions with the concentration 
inequality based bounds in the context of compressed sensing. 
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Fig. 1. Distribution of the smallest eigenvalue for Wishart real matrices VV’s(m, I), m = 400, s = 10. Comparison between 
the exact (dashed line), the gamma approximation (61) (dotted), and the concentration inequality bound (75) (solid). 



Fig. 2. Distribution of the largest eigenvalue for Wishart real matrices >Vs(m, I), m = 400, s = 10. Comparison between the 
exact (dashed line), the gamma approximation (60) (dotted), and the concentration inequality hound (74) (solid). 
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